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Abstract 

In this paper we consider Runge-Kutta discontinuous Galerkin (RKDG) schemes 
for Vlasov-Poisson systems that model collisionless plasmas. One-dimensional systems 
are emphasized. The RKDG method, originally devised to solve conservation laws, is 
seen to have excellent conservation properties, be readily designed for arbitrary order of 
accuracy, and capable of being used with a positivity-preserving limiter that guarantees 
positivity of the distribution functions. The RKDG solver for the Vlasov equation is 
the main focus, while the electric field is obtained through the classical representation 
by Green's function for the Poisson equation. A rigorous study of recurrence of the 
DG methods is presented by Fourier analysis, and the impact of different polynomial 
spaces and the positivity-preserving limiters on the quality of the solutions is ascertained. 
Several benchmark test problems, such as Landau damping, two-stream instability and 
the KEEN (Kinetic Electro static Electron Nonlinear) wave, are given. 

Keywords: Vlasov-Poisson, discontinuous Galerkin methods, recurrence, positivity-preserving, 
BGK mode, KEEN wave. 



1 Introduction 

The Vlasov-Poisson (VP) system is an important equation for modeling collisionless plasmas, 
one that possesses computational difficulties of more complete kinetic theories. Thus, it serves 
as an important test bed for algorithm development. The VP system describes the evolution 
of / = f(x, v, t), the probability distribution function (pdf) for finding an electron (at position 
x with velocity v at time t) with a uniform background of fixed ions under a self-consistent 
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electrostatic field. In particular, the non-dimensionalized VP system (with time scaled by the 
inverse plasma frequency u p 1 and length scaled by the Debye length Ad) is given by 



d t f + vV x f-E-V v f = Ox(0,T] 
-A s $ = l- f fdv n x x (0,T] 



E = -V x $ Q x x(0,T}. 

Here the domain Q = Q x x IR n , where Q x can be either a finite domain or IR n . The boundary 
conditions for the above systems are summarized as follows: / — > as \x\ — > oo or |i>| — > 
oo. If Q x is finite, then we can impose either inflow boundary conditions with / = f m on 
Tj = {(x, v)\v ■ v x < 0}, where v x is the outward normal vector, or more simply impose 
periodic boundary conditions. For simplicity of discussion, in this paper, we will always 
assume periodicity in x. Also, we add that when the VP system is applied to plasmas the 
total charge neutrality condition, J n (J Rn f dv — l) dx = 0, is imposed. 

The following physical quantities associated with this system are related to its conservation 
properties: 



charge density p(x,t) = / f(x,v,t)dv, 

momentum density j(x,t) = / vf(x,v,t)dv, (2) 

kinetic energy density ^(x,t) = - f \v\ 2 f(x, v , t) dv , 

electrostatic energy density !; e (x,t) = -\E(x,t)\ 2 . 

Indeed, it is well-known that the VP system conserves the total electron charge j Q p(x) dx, 
momentum J n j(x)dx, and energy j Q (^fc(x) +£ e (x))dx. Moreover, any functional of the 
form J n G(f) dxdv is a constant of motion. In particular, this includes the k-th order invariant 
Ik = J n f k dxdv and the entropy S = — f n fln(f) dxdv. Sometimes the functional I2 is also 
called the enstrophy, and all of these invariants are called Casimir invariants (see, e.g., [34J). 

The VP system has been studied extensively for the simulation of collisionless plasmas. 
Popular numerical approaches include Particle-In-Cell (PIC) methods [SI [21] , Lagrangian par- 
ticle methods [U [TTj , semi-Lagrangian methods [HI 05], the WENO method coupled with 
Fourier collocation [5H], finite volume (flux balance) methods [TBI EE] , Fourier- Fourier spec- 
tral methods [271 128] . continuous finite element methods [HJ [50], among many others. In this 
paper, we will focus on the discontinuous Galerkin (DG) method to solve the VP system. The 
original DG method was introduced by Reed and Hill |42j for neutron transport. Lesaint and 
Raviart [32] performed the first convergence study for the original DG method. Cockburn and 
Shu in a series of papers [HI [13J [121 El US] developed the Runge-Kutta DG (RKDG) method 
for hyperbolic equations. The RKDG methods have been used to simulate the VP system in 
plasmas by Heath, Gamba, Morrison and Michler [231 [22] and for the gravitational infinite 
homogeneous stellar system by Cheng and Gamba [H] . Theoretical aspects about stability, ac- 
curacy and conservation of those methods are discussed in [221 [23] and more recently in [31 [2] 
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for energy conserving schemes. Such methods have excellent conservation properties, can be 
readily designed for arbitrary order of accuracy, and have the potential for implementation on 
unstructured meshes with hp— adaptivity. To ensure the positivity of the solution, one can use 
a maximum-principle-satisfying limiter that has been recently proposed by Zhang and Shu in 
[52] for conservation laws on cartesian meshes, and later extended on triangular meshes |56j . 
This limiter has been used to develop positivity-preserving schemes for compressible Euler 
[531 [55], shallow water equations [H], and Vlasov-Boltzmann transport equations [10] . It has 
also been employed recently in the framework of semi-Lagrangian DG methods (131 [H] for the 
VP system. 

The scope of the present paper is as follows: we focus on a detailed study of the RKDG 
scheme for the Vlasov equation from both the numerical and analytical points of view. Since 
we are only considering one-dimensional problems, we use the classical representation of the 
solution by Green's function to compute the Poisson equation; therefore, the electric field is 
explicitly given as a function of the numerical density. This removes all discretization error 
from the Poisson equation and lets us more accurately investigate our DG solver for the Vlasov 
equations. We rigorously study recurrence, which is an important numerical phenomenon that 
commonly appears with many solvers. We use Fourier analysis and obtain eigenvalues of the 
amplification matrix, and then investigate the impact of different polynomial spaces on the 
quality of the solution by examining conserved quantities as well as convergence to BGK states 
for some choices of initial states. We consider benchmark test problems such as simulations 
of Landau damping phenomena for the linearized and nonlinear Vlasov Poisson systems, two- 
stream instability, and their long time BGK states and the formation of KEEN waves, both 
for the nonlinear system as well. 

The remaining part of the paper is organized as follows: in Section [2j we describe the 
numerical algorithm and summarize its conservation properties. In Section [3j we study the 
recurrence phenomena that occurs for linear Vlasov type transport equations discretized by 
DG methods with various polynomial orders. Sections 4A_ and <L2_ are devoted to discussions of 
simulation results for the linearized and nonlinear VP system, respectively, for diverse choices 
of initial data and external drive potentials. We conclude with a few remarks in Section [5] 



2 Numerical methods 

In this section we first describe the proposed DG numerical algorithm and then discuss some 
of its basic conservation properties related to the quantities of ([2]). This is done for both the 
fully nonlinear VP system of ([I]) and the linearized VP system obtained by linearizing about 
the homogeneous equilibrium f eq {v), with corresponding vanishing electric equilibrium field. 
Periodic boundary conditions in x are assumed. 

Thus, setting f(x,v,t) = f eq (v) + 5f(x,v,t) and expanding the system to first order ap- 
proximation, the perturbation 5f satisfies the Linear Vlasov-Poisson (LVP) system, 

dtf + v- V x f = E- V v f eq Qx(0,T} 

A x $= f fdv n x x(0,T] (3) 
E = -V x <$> n x x(0,T], 
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where 5f has been replaced by / to ease the notation. We find it convenient and efficient to 
intertwine the discussion of our algorithms for the VP and LVP systems. To avoid confusion 



in Section 2.1 we underline the words linear and nonlinear, signaling where discussions specific 



to each apply. 

2.1 Numerical algorithm 

For one-dimensional problems, we use a mesh that is a tensor product of grids in the x and 
v directions, because this simplifies the definitions of the mesh and polynomial space for the 
Poisson equation. Specifically, the domain Q is partitioned as follows: 

= xi_ <X3 < ... < x Nx+ i — L, -V c = vi < va < . . . < v Nv+ i_ = V c , 

where V c is chosen appropriately large to guarantee f(x,v,t) = for \v\ > V c . This is a 
reasonable assumption, because of the well-posedness of the one-dimensional Vlasov-Poisson 
system as indicated in [21 J. The grid is defined as 



= X [Vj_l,V j+ i], 

Ji = [x i - 1 /2,x i+1/2 ], Kj = [v j - 1 /2,v j+1 / 2 ], i = l,...N x , j = l,...N, 



where Xi = h(x { _i + x i+ i) and Vj = Mvj^u 2 + fj+1/2) are center points of the cells. 

We will make use of several approximation spaces. For the x-domain, we consider the 
piecewise polynomial space of functions £ : Q x — > K, 

Z l h = {V- eke^(Ji). i = l,...N m }, 

where P l {Ji) is the space of polynomials in one dimension of degree up to I. For the (x,v) 
space, we consider two approximation spaces of functions <fi, <p : Q — > R, 

Vk = {0 = e Q?(I id ), i = l,...N x , j = 1, . . . N v } 

and 

W l h = W: <p\j iti e P'(/ M -), i = l,...N x , j = 1, . . . N v } , 

where Q l (I itj ) = P\Ji) ® P l (Kj) = span{x h v l \ VO < h < /, < l 2 < 1} denotes all 
polynomials of degree at most I in x and v on I id , and P z (/ij) = span{x l1 v l2 ,W0 < l\ + 1 2 < 
I, li > 0, l 2 > 0}. These two spaces are widely considered in the DG literature for multi- 
dimensional problems. A simple calculation shows that the number of degrees of freedom of 
Q is (I + l) 2 . For I > 1 this is larger than the number of degrees of freedom of P z (ijj), 
which is (1 + l)(/ + 2)/2. 

First we describe the RKDG scheme for the linear Vlasov equation. We seek fh(x, t) E V} t 
(or W l h ), such that 



(fhjtVhdxdv - / vf h (<p h ) x dxdv + / (v f h tp h v dv 

Jl itj .lis, 2 ' 

-/ (Vfh<pt)i-i,v dv = / E h f' eq (p h dxdv (4) 

JKi Jin 
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holds for any test function ifh{x,t) G V h (or W l h ). Here and below, we use the following 
notations: Eh is the discrete electric field, which is to be computed from Poisson's equation, 

= lim e-x) <Ph(x i+1 / 2 ± e,v), {Vh)% j+1 / 2 = lim^o <Ph(x, v j+1 / 2 ± e), and vfh is a nu- 
merical flux. We can assume that in each Kj, v has a single sign by properly partitioning the 
mesh. Then, the upwind flux is defined as 

~ I vfh if v < in Kj • 1 j 

The scheme for the nonlinear Vlasov equation is similar. We seek fh(x,t) G V h l (or W l h ), 
such that 



{fh)t<fh dxdv - / vf h (ip h ) x dxdv + / v dv- „ 

+ / E h f h (tp h ) v dxdv - / (E h f h (p-) j+ idx+ / (E h f h ip+) j_i dx = (6) 

holds for any test function ip h {x,t) G (or The upwind flux for vfh has been defined 

in ([5]) and the new flux needed for the nonlinear case is given by 

hU ~ l if Ij t Ehdx > • (7j 

The above descriptions coupled with a suitable time discretization scheme, such as the 
TVD Runge-Kutta method jS], completes the RKDG methods. For example, suppose the 
semi-discrete schemes in Q and (|6| are written in the compact form 



(h)t<Ph dxdv = Uijifh, E h , <Ph) 

i 

where for the linear Vlasov of Q 

K l i,j(fh,E h ,(Ph) = I vf h (iph) x dxdv - I {vf h (fh) i+ i v dv 



+ / (vfh<Ph)ir-Lv dv + / E h f' eg (p h dxdv, 
Jk 3 Juj 

while for the nonlinear Vlasov of (J6J) 
Hlf in (f h ,E h ,cp h ) = I vf h (<f h ) x dxdv- I (ff h <Ph) i+hv dv+ [ (vfhVttt-frdv 

*i 3 3 ^3 

- I E h fh(Vh)vdxdv + / (E h fhVh)x,j+\dx - \ (E h f h ipl) xJ _i dx . 

J 14 A J Ji J Ji. 
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The third order TVD-RK method implements the following procedure for going from t n to 

t n+l. 



f f®<p h dxdv = lf fttp h dxdv + i f fV><p h dxdv + E h\ <Ph) , (8) 

( fl+Vfc dxdv = I / ft<p h dxdv + | / /<* V fc ^ + ^^(/f , 4 2) , p h ) . 

Poisson's equation is used to obtain J5£, -E^, and E h 2 \ Beyond periodicity, we need to 
enforce some additional conditions to uniquely determine For example, we can set one end 
of the spatial domain to ground, i.e. set $(0, t) = 0. In the one-dimensional case, then the 
exact solution can be obtained if we enforce $(0) = <&{L). For the nonlinear system we obtain 

®h = j j Ph(z, t) dzds - — - C E x, 



(9) 



o jo 



where Ce = —L/2 + J Q L J* ph(z, t) dzds/L, and 

E h = -($ h ) x = C E + x - / p h {z,t)dz, 

Jo 

while for the linear system Poisson's equation gives 

®h= / Ph(z,t)dzds - C E x, 
Jo Jo 

where C E = J Q L J Q S Ph(z,t) dzds/L, and 

E h = -($ h ) t = C E - [ p h (z,t)dz. (10) 



o 



From (9) and (10), we see that if f h G V{ (or W l h ), then p h = J^ v f h dv G Z l h ; hence, E h G 
Z l h +1 f\C°. The above approach uses the classical representation of the solution by Green's 
function and will be referred to as the "exact" Poisson solver. It is valid only for the one- 
dimensional case. For higher dimensions, a suitable elliptic solver needs to be implemented, 
such as those discussed in [23]. Here we use the exact solver to entirely eliminate discretization 
error from Poisson's equation and, thereby, spotlight the performance of the Vlasov solver. 

Below we describe positivity-preserving limiters, as summarized in [31]. We only use such 
a limiter to enforce the positivity of fh for the nonlinear VP system. For each of the forward 
Euler steps of the Runge-Kutta time discretization, the following procedures are performed: 

• On each cell Iij, we evaluate T^j := min^)^ ^ fh{x,v), where S it j = (Sf <8> <SJ) [j{Sf <g> 

Sj), and Sf,Sj denote the (I + 1) Gauss quadrature points, while Sf,Sj denote the 
(I + 1) Gauss-Lobatto quadrature points. 
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• We compute f h (x, v) = 9(J h (x, v) - (fhUj) + (fh)i,j, where (fh)i,j is the cell average of f h 
on Jjj, and 9 = min{l, \(fh)i,j\/\Tij — (fh)i,j\}- This limiter has the effect of maintaining 
the cell average, while "squeezing" the function to be positive at all points in S it j. 

• Finally, we use fh instead of fh to compute the Euler forward step. 
This completes the description of the numerical algorithm. 



2.2 Scheme Conservation properties 

In the following, we will briefly review and discuss some of the conservation properties of the 
above RKDG scheme for the nonlinear VP equations without the positivity-preserving limiter. 
Some of those results have been reported in [22] and [3]. 

Proposition 1 (charge conservation) For both the and W l h spaces, 

Y^HZf in (f h ,E h ,l) = e(f h ,E h ,l) 

which implies 

W f^dxdv = W ftdxdv+^At(e(fl 2 \E% l \l) 

+\®(fi 1 \E£\i) + 1 l e(fZ,Eii)) 

for the fully discrete scheme Q). Here, 

e(/ h , E h , <p h )=^2 / ( E hfhVh) x , Nv+ i dx-^ ( E hfhVh) x ,i 2 dx 

i * i i 

denotes the contribution from the phase space boundaries located at v = ±V C , and should be 
negligible if V c is chosen large enough. 

Remark: Charge conservation (or mass conservation or probability normalization as it 
is sometimes called) states that the total charge will be preserved on the discrete level up to 
approximation errors associated with the phase space boundaries. The proof is straightforward 
and, therefore, omitted. The same conclusion can be proven for the linearized system. The 
positivity preserving limiter does not destroy this property because it keeps the cell averages 
unchanged. 

Proposition 2 (Semi- discrete L 2 stability - enstrophy decay) For both the V[ and W l h spaces, 
n^ih, E h , fh) < 0. Hence, 
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The proof, for an arbitrary field Eh, can be found in [ID] . Theorem 4, which applies directly 
here by setting the collisional form Q a = in that proof. 

For the remainder of this subsection we will assume the DG solution satisfies the ve- 
locity boundary conditions fh(x,±V c ,t) = 0. This is a reasonable assumption when V c is 
large enough. In particular, this will guarantee exact charge conservation, which implies that 
j ph{x,t)dx is constant in time t. Therefore, using the definition of Eh in (9), we can obtain 
periodicity in Eh, i.e, £^(0) = Eh(L). Without this assumption the propositions below contain 
multiple boundary terms and the proof becomes technical. 

Proposition 3 (Momentum conservation) Assuming fh(x,±V c ,t) = 0, for both the and 
W l h spaces when I > 1, ^\ ■ H-j nhn (A, Eh, v) = 0, which implies 



J2 f f£ +1 vdxdv = J2 [ Kvdxdv 



for the fully discrete scheme. 



vf h v)i_i v dv 

Ki 



Proof. Choosing cph — v in ([6j, we have 

J2n?f in (f h ,E h ,v) = V / vf h (v) x dxdv- [ (ff h v) i+ i dv+ f 

- I E h fhdxdv+ / (E h f h v) X)j+ i dx - / (E h fhv) x>j _idx 
= - ^ / E h f h dxdv = - ^2 / E hPh dx , 

and using the exact Poisson solver together with the periodicity of Eh and $^ yields the 
following: 

^2 / E hphdx = / ^(pfe - 1) dx + / £ fo da; 

= -J2 / + / Sfcdar 

= -(El(L)-El(0))/2-^L) + m = 0, 

which completes the proof. □ 

Remark: The above proof holds for the linearized system as well. Note, however, it 
relies on the use of the exact Poisson solver. For a full numerical DG Poisson solver, such as 
that developed in [23] for the discretization Poisson equation, exact momentum conservation 
remains true, as was proven in [22] by means of the DFUG method developed there for dealing 
with the discretized Poisson equation. However, the positivity-preserving limiter we use here 
will destroy this property because it was not designed to conserve the numerical momentum. 
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Proposition 4 (Semi- discrete total energy equality) Assuming fh{x,±V c ,t) = 0, for both the 
V h l and W l h spaces when I > 2, 



d 

It \ 2 



\Y, \ dxdv + \Y,I e ZW dx ) = Mfh, *h) = Mh - f, *h - P(*h)) 



where the operator A(w,u) := [ wu x v ~ (w)tU) dxdv, and P is any projection such 

that P($ h ) e Z[ and P{$ h ) = $ h at x i+1/2 , for i = 0,1, ... , N x . 



Proof. Choosing t^ h = v 2 /2 in (|6| yields 

y^l, J fav 2 dxdv + ^2 I E hfhV dxdv = 



d 
dt 



and 



d 
dt 



Y^\l E 2 h {x)dx = [ E h (E h ) t dx = -J2 I {®h) x {Eh)tdx 

i J Ji i J Ji i J Ji 

= Yl / ®h(E h ) xt dx = Y / $hO--Ph)tdx 

= ®h{ph)idx = -^2 $h(fh)tdxdv, 

where in the second line, we have used the periodicity and continuity of Eh and Therefore, 
we have proven that 



d 
dt 



f h v 2 dxdv+±J2 I E 2 h (x)dx) =A(f h ,$ f . 



On the other hand, upon choosing iph = P(&h) in ^ and using the periodicity and continuity 
of P($ fc ), we can verify that A(f h , P($ /l )) = 0. The exact solution / obviously satisfies 
A(f, §h — P($/j)) = from the continuity and periodicity of $^ — P($^), and therefore we 
are done. □ 

The above proof indicates that the variation in the total energy will be related to the 
error of the solution, fh — f, together with the projection error, §h — P(&h)- I n [221 E3], 
error estimates for DG schemes with NIPG methods for the Poisson equations are provided 
for multiple dimensions. In [3], optimal accuracy of order I + 1 for the semi-discrete scheme 
with Q 1 spaces has been proven under certain regularity assumptions. We remark that in 
[3] conservation of the total numerical energy is proven when the Poisson equation is solved 
by a local DG method with a special flux. Unfortunately, no numerical simulations of linear 
Landau damping or of the nonlinear VP system, such as those done in [23] or in Section [4] of 
this present manuscript, have been performed up to this date by the scheme proposed in [3], 
so a comparison is not possible. 
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3 On recurrence 



In this section we study recurrence, a numerical phenomenon that is known to occur in simu- 
lations of Vlasov-like equations. Its study is important because it provides information about 
the numerical accuracy of a scheme. Recurrence was observed in the semi-Lagrangian scheme 
of Cheng and Knorr [8], where a simple argument for its occurrence was provided. In this 
section, we carry out a detailed study of recurrence for the DG method. 

We study recurrence of our algorithm applied to the linear advection equation ft + vf x = 
on [0, L = 2ir/k] x [— V c , V c ], since it is tractable and contains the basic recurrence mechanism; 
results for the full Vlasov system tend to be quite similar. The initial condition we consider 
is fo(x,v) = A cos(kx) f eq (v) , and the equilibrium distribution f eq iy) is taken to be either the 
Maxwellian or Lorentzian distribution, viz. 

fM = 4^e~^ or / ' 1 



2vr irv 2 + l 

For the Maxwellian equilibrium, J'm, we take V c = 5, and for the Lorentzian equilibrium, /x, 
we take V c = 30. 

The exact solution for the advection equation is f(t, x, v) = fo(x — vt, v). Hence, a simple 
calculation shows p(x,t) = Acos(kx)e~ k * I 2 for the Maxwellian distribution; similarly, for the 
Lorentzian, p(x,t) = Acos(kx)e~ kt . Thus, we see how the density for each case should decay 
to zero. The failure of decay and the occurrence of recurrence noted for the semi-Lagrangian 
scheme of [8] stems from the finite resolution in the velocity space and, indeed, the recurrence 
time depends on the mesh size in v. 

To be specific, we repeat the definition of DG scheme for this equation, which amounts to 
(p| with E h set to zero: we find fh{x,t) £ V£ (or , such that 



(fh)tfhdxdv- vf h (ip h ) x dxdv+ J (vf h tp h ) i+ i v dv- j (vf h ip^) i _i v dv = (11) 



holds for any test function (p h (x,t) £ V* (or W k ). Again vfh is the upwind numerical flux 
of ([5]). In the analysis below, we always assume time to be continuous, because recurrence is 
mainly a phenomenon that comes from the spatial and velocity discretization. 



3.1 The case of / = 

For the piecewise constant case, the DG method is equivalent to a simple first order finite 
volume scheme and we can derive rigorously the behavior for p. Suppose we define fh = fij on 
cell Iij, and assume uniform grids in both directions. Moreover, we assume N v to be even for 
simplicity. With this assumption, the location of the cell center is Vj = (j — jV " 2 +1 ) Av . Now 



(11) simply becomes 



m. + v ti+ij-te = if u, < 0. 

at J Ax J 



(12) 
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The initial condition chosen is clearly equivalent to fij(0) = Re (Ae lkXi f eq (vj)) . We prove that 
the scheme above gives 

f lJ (t) = Re(Ae ikx ^ t f eg (v J )) (13) 



where Sj is given in (14) below 



Upon plugging (13) into (12), we have 



f I 1— e~ ikl 



-fx 



if Vj > 



Hence, 



— v 



1 _ e -ikAx _ cos(fcAx)-l 



J Ax 

gifc Ax 



J 3 Ax 



J Ai J — 



J Aj 

which can be summarized as 



cos(fcAx) — 1 s'm(kAx) • -r . r\ 

-Vj — S; — Vj — 7- — -1 it Vj < U , 

J Ax J Ax J ' 



Sj = \v j{ 



.cos(kAx) — 1 sin(fcAx) . 



Ax 



Ax 



1 . 



(14) 



Therefore, the real part of Sj is always negative, this means the magnitude of will always 
damp, but because of the j-dependence it does so at different rates for different cells. Since 
the density 



P( x i) = ^2 k Av = Re 



^ j^ e ikxi+Sjt 



feq(Vj) AV 



the density will damp at a rate between 



Av cos(kAx) — ! 



and 



V c —Av cos(kAx) — 1 



Another impor- 



2 Ax 2 Ax 

tant fact is that recurring local maxima of the density will have a period Tr that satisfies 



Av sin(fcAx) 



ir = 71 . If we define k' = kia ^- L > then Tr 

2 Ax n Ax ' n 

coincides with the recurrence time obtained in [8]. 

Next we compare the above theoretical prediction with numerical results. In all of the 
calculations below, we take A — 1, k — 0.5 and the mesh size to be 40 x 40. In Figure [TJ we 
display results for numerical runs using piecewise constant polynomials and time discretization 
using TVD-RK3. (We use the third order method to minimize the time discretization error.) 
We plot p max (t) = m&x x p(x, t) in Figure [TJ First, we notice the pattern of p max has the 
expected periodic structure with damping for both Maxwellian and Lorentzian equilibria. 
For the Maxwellian distribution, a simple calculation yields Tr = 50.47. Similarly, with the 
formulas above, the smallest damping rate is —0.49 x 10~ 2 , while the biggest is —9.3 x 10~ 2 . 
For the Lorentzian distribution, Tr = 8.41, and the smallest damping rate is —2.94 x 10~ 2 , 
while the biggest is —5.58 x 10 _1 . From Figure [TJ by using the second to the fourth peak, the 
actual computed value of Tr for Maxwellian is 50.32 and the damping rate is —1.02 x 10~ 2 ; 
while for Lorentzian, from the second to the tenth peak, Tr is 8.40 and the computed damping 
rate is— 3.19 x 10~ 2 . These numbers agree well with the theoretical prediction. 

In conclusion, our analysis explains the behavior of the first order DG solution. At t — Tr, 
the numerical density obtains a local maximum; hence, clearly at this time the numerical 



s'm(kAx) 



2tt 
k'Av ' 



When Ax — > 0, k' — > k, and this 
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solution can no longer be trusted. The numerical decay deviates from the theoretical decay 
well before t = Tr. To achieve a larger Tr, according to the formula, we can refine Av. On 
the other hand, refining Ax will not change Tr by much. 




n1(T 
E 

I 

CL 

1 ' 



100 

t 



fa) Maxwellian, P° 



(b) Lorentzian, 



Figure 1: Computations of the advection equation for piecewise constant polynomials showing 
local maxima of the density p max as a function of time. The mesh is 40 x 40 with Ax = n/10. 
For the Maxwellian equilibrium Av = 1/4, while for the Lorentzian equilibrium Av = 3/2. 



Remark: Using the same methodology, it is easy to perform a similar analysis for any type 
of finite difference method. The real part of Sj will be negative if there is numerical dissipation, 
and the imaginary part will always approximate Vjk due to the differential operator v-^. This 
means that for such schemes, the recurrence time Tr will always be close to j^-. 



3.2 Higher order polynomials 

In this subsection, we consider higher order polynomials. For the space, it takes four point 
values in each cell to represent a Q 1 polynomial. This technique was developed in [FT] for 
analyzing piecewise linear DG solutions in one dimension. As in Section 3.1 we use a uniform 
mesh, 

of v > only, then vfh 



Axi = Ax and Avj = Av. 



Without loss of generality, we consider (JTTj) for the case 

+ 1. 



vf h , which means we only consider cells ? with j > 



Nv 
2 



In each computational cell Iij, we can always use the following form to represent fh- 



fh = fi-y+i Xi(x,v) + fi-y-i X2{x,v) + X3(x,v) + /, 



i+hj- 



Xi(x,v), 
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where 



Xi(x, 


v) 


= -4 


X2(X. 


v) 


-«( 


X3(X, 


v) 


-«( 


Xa(x. 


v) 


= -4 



1 

Axi ~ 4 

X — X,; 1 



V — Va 1 
+ - 

Av< 4 



Ax,-, A' 
Ax,- + 4 



+ - 

Axi 4 



V-Vj 1 

4 

1 



Av, 



Aw,- + 4 



are the basis functions in Q 1 and /i±i/4,j±i/4 = fh( x i±i/4i v j±i/4) are the point values. By 
choosing the test function in (11) to be ifh — Xe.t ^ — 1)2,3,4, we obtain four relations. 
Letting f tj = (/i_i/ 4 j+i/4, /i-i/4,j-i/4, /i+i/4j+i/4, /i+i/4,j-i/4) T , then corresponding to each of 
the four terms in (11), we have 



eft 



Bfij + Cfi j - Df, 



0. 



where 



M 



AxAv 



144 



/49 -7-7 1 \ 

-7 49 1 -7 

-7 1 49 -7 

V 1 -7 -7 49/ 



B 



C 



Av 
~Y2 



Av 
48" 



(-{2Av + 7^-; 



and 



V 3 

2Av + 7fj 

V ~ v i 
( 2Av + 7vj 

-6 Aw - 21^- 

3Vn 



/-6Av - 21vj 



2Av - 7vj 



~ V 3 

-2Av + 7vj 

-2Av + 7vj 

3vj 
6Av - 2lvj 



■(2Av + 7vj 



2Av + 7v; 



-(6Av + 21vj 
3vj 

18 Av + Q3vj 
-9 Vj 



D 



Av 
48" 



3vj 
2Av + 7vi 



3vj 
6Av - 21vj 

\ -vj -2Av + 7vj 

After simple algebraic manipulation, we obtain 

dfjA Av 



18 Av + 63vj 

-9vj 
-6Av - 21vj 
3vi 



Vj \ 

2 Av - 7vj 

-2Av + 7vjJ 

3vj \ 
QAv - 21vj 
-9vj 
-18Av + 63vjJ 



-9vj 
-18Av + 63uj 
3vj 

6AV-21VJ- ) 



dt ~ Ax ( S ™fv +T ™f*-u) 
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where 



/ 49 _ 
' 96 



7m 



Sm. — 



7 



96 

49 7m 



32 



3m 



32 
7 3r 



T 



and m = 2j — N v — 1 
matrix is given by 



/- 

V 

1,3,5 



96 

77 , 11m 
96+ 8 

96 

35 5m 



96 



11 



77 96 
77 i m 



1_ 

96 



_5_ 
96 



96 



96 

35 5m 



32 

21 9m 

32 3 8 

_ 32 
15m 



35 
32 



+ 



96 



1 



„ 32 

7 3m 



32 8 
_3_ 

32 

21 9m 
32 8 
_5_ 
32 

35 I 15m. 

32 8 



96 



7 96 

]_ \ rn 
96 8 



- 32 
7 3m 



32 8 
~~ 32 32 8 

are positive and odd integers. Therefore, the amplification 



/ 



Gi 



Av 



(S m + T m e~ lkAx ). 



3 Ax 

With the initial condition /y(0) = Re(Ae lhXi T), where 



T 



^-ikAx/A 



feq(v j+ l),e 



—ikAx 



ikAx 



AkA 



x/A j 



it is clear that the general expression for the numerical solution is 



f ij (t)=Rele ikXi Y t a a V a e 



Vat 



Here r] a are the eigenvectors of Gj with V a the corresponding eigenvectors, a a are constants 
such that fij(0) = T, and all these quantities are dependent on j (or equivalently m = 
2j — N v — 1). The collective behaviors of the eigenvalues r] a will influence the behavior of the 



density as a function of time. We focus on the matrix A r 
algebraic manipulation can be seen to have the form 

Am = W (g) V, 



S m -\- T m c 



-ikAx 



, which with some 



where W and V are the following 2x2 matrices: 



W 



V = 



—3m — \ 
_ i 

4 

2 ^ 24 l 



1 

4 

—3m 



2 24 * 



2 Ij 
2 + 8 2 



and i = e lfeA:E — 1 = —ikAx + 0(Aa; 2 ) 
formulations of the mesh and the space 
+ V9 + 12i + i 2 )/6 



obtaining A] 



This nice structure is due to the tensor product 
l . We compute the eigenvalues of the matrix V, 
UkAx - ±k 2 Ax 2 + 0(Ax 3 ) and A 2 = (3 + % + 



6 12 

0(Ax 2 ), and the eigenvalues of W, obtaining —3m + \/3. 



^9 + 12i + i 2 )/6 = 1 - \ikAx 
Hence, by simple linear algebra, the four eigenvalues of the matrix A m are obtained 

/6 = (-3m - V3)X 2 \ 
& = (-3m + v / 3)A 2 
e 3 = (-3m- v / 3)A 1 
\^4 = (-3m + y/3)\J 



14 



It is easy to show that the eigenvectors corresponding to these eigenvalues are independent of 
m, since the eigenvectors of V and W are independent of m. We conclude that the eigenvalues 



of Gj are 



'-3m - V3)\ 2 Av/Ax\ 
[-3m + V3)X 2 Av/Ax 
'-3m - v^AiAw/Ax 
'-3m + ^XiAv/AxJ 



and therefore, 



^2a a V a e^ 



-3m(X 2 Av/Ax)t 



aiVi e 



-V3(X 2 Av/Ax)t 



+ a 2 V 2 e 



V3(X 2 Av/Ax)t 



a=l 



+e 



-3m(XiAv/Ax)t ^y^ ^^(X^v/Ax^ + a ^ ( ,V3(X 1 Av/Ax)t^ 



Since rji and rj 2 have nontrivial negative real parts, the damping for those two modes will be 
strong. Consequently, the main behavior of the density will be dominated by the eigenmodes 
of rj 3 and 7? 4 . Recall 

P (x i± i, t) = ^2J fh(x i± i, v, t)dv = J2(fi±y+\ + fi±\ Av ' 



and, therefore, the behavior of p(xi±i/4,t) is dominated by 



-3m(X 1 Av/Ax)t 



r > ( , V3{\!Av/Ax)t _|_ ^ ^(XiAv/Ax)^ 



where C3 and C4 are constants that do not depend on m. Since Ai = \ikAx—^k 2 Ax 2 +0{Ax 



we have — 3m(AiAt>/Ax 
to that of Section 
of all modes will 



kAv „ 



-mi- 



mAvAxk 2 



3.1 



2 """ 4 
for the piecewise constant case, when t 



\-0(AvAx 2 ). Hence, with an argument similar 



Tr = t^ 1 the imaginary parts 
return to mni, and this will correspond to a local maximum of p ma x as a 
function of time. The remaining term, ce~' v3tkAv/m ^ ^ e V3ikAv/6t ^ correS p 0nc l s to the envelope 
of the wave, and the negative real part of the eigenvalues indicates numerical dissipation. 

In Figure |2l we plot the evolution of p m ax as a function of time for the Q 1 and Q 2 spaces. 
From Figures ^[a) and |2](b), we observe the behavior predicted by our analysis for the Q 1 
elements. From Figures^2](c) and[2|d), we find that the solutions using the Q 2 polynomials 
share similar structures, except that small oscillations can be observed for the Maxwellian 
case. Also we note that the Q 2 discretizations can follow the exact solutions longer in time, 
in the sense that the minimum value achieved before p max starts to deviates from the exact 
solution is on the order of 10 -6 compared to 10~ 4 in the Q 1 case. This is expected due to the 
higher order accuracy of the scheme. For the Q 2 polynomials, we deduce that the amplification 
matrix G is a 9 x 9 matrix. Thus, for this case there, there will b e nine eigenvalues and more 
modes than for the Q 1 space. In Table [TJ we verify the recurrence time Tr numerically; good 
agreement between the predicted values and the observed values are seen. 

Note, the trace Tr(S m ) = -4m,m = 1,3,5. . ., while a similar calculation for cells when 
v < yields Tr(S m ) = 4m, m = 1,3,5.... Therefore, we conclude that our semi-discrete 
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Table 1: The location of local maxima of the density p max compared with the predicted 
recurrence time Tr. The Tr values for the Maxwellian equilibrium are computed using the 
average of the first three peaks, while for the Lorentzian they are computed using the average 
of the first seven peaks. 





Predicted T R = 


Numerical value of Q 1 


Numerical value of Q 2 


Maxwellian 


50.26548245743669 


50.265482457450 


50.265482457450 


Lorentzian 


8.37758040957278 


8.37787960887962 


8.37787960887760 



algorithm has an incompressible vector field and thus possesses a version of Liouville's theorem 
on conservation of phase space volume. Liouville's theorem for finite difference and Fourier 
discretization of fluid and plasma equations is well known and has been used in statistical 
theories of turbulence (see e.g. [401 ED ESI ES])- We also note that we have performed the 
analysis for the semi-discrete DG schemes. For fully discrete RKDG schemes, one could 
use a similar method, as proposed in [57] for the wave equation, to write the fully discrete 
amplification matrices, but we do not pursue this in this paper. 

Finally, we note in closing this section that a similar analysis using the P 1 elements yields 
a 3 x 3 matrix; however, this basis does not yield the nice form possessed by Q 1 because of 
the loss of the tensor structure. Figure [3] shows the temporal behavior of p max using the P' 
elements. Observe that, although the local maxima still are located near Tr = there 
appear to be several small local maxima instead of one main maximum, and overall the long 
time dissipation seems to be stronger than that for Q 1 . We also noted that the P 2 basis follow 
the exact solution longer than P , but shorter than Q 2 cases, because for P 2 elements, the 
minimum value that the solution achieves before it deviates from the exact solution is on 
the order of 10 -4 . In summary, we conclude that increasing the polynomial order does not 
change Tr by much. However, higher order accuracy seems to improve the time that the 
numerical solution can follow the exact solution. For Q 1 elements, the amplification matrix 
can be written as a tensor product of two small matrices, and this made possible our direct 
analysis for the recurrence time. For P z elements, we lost the tensor structure, and the solution 
is more dissipative. 

4 Vlasov numerical results 

Now we turn to some numerical tests of our method for both the VP and LVP systems. 
For the LVP system we consider the standard tests of linear and nonlinear Landau damping, 
which have been studied in many references in the contexts of various numerical techniques 
since [S] (see [23] for an extended list), but we also consider a test that heretofore does not 
appear to have been done, viz., we monitor the linearized energy that is conserved by the LVP 
system [30] EHl [35]. Similarly, for the nonlinear VP system we consider the standard tests 
of nonlinear Landau damping and a symmetric version of the two-stream instability (also see 
[23] for references). In addition, for the VP system we consider an example that is initialized 
by a driving electric field, resulting in a dynamically accessible initial condition as described 
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in [36], EZl [38] , which has been observed to approach nonli near BGK [5] states that have been 
termed KEEN waves in Ref. [H [25] (see also [22]ll6]). 

4.1 Linearized VP system 

Associated with the LVP system of ^ is the well-known plasma dispersion function |20j . 

1 f + °° /» 



e(k,u J ) = l-- ^jrdv, (15) 
k 2 J-oo v-u/k 

which (with the appropriate choice of contour) will be used to benchmark the accuracy of 
the Landau damping rate and oscillation frequency obtained from our DG solver with choices 
for the various polynomial spaces. The LVP system conserves not only the total charge and 
momentum, but also the linear energy J3U1 EH1 ES] , which is defined as 

i r v f 2 if 

H L = -- I ^L-dxdv + - I E 2 dx. (16) 



2 Jfl f'eq 2 

As noted above, we monitor this quantity and check for its conservation. In addition, we 



monitor the shift of energy to the first term of ( 16 ) as the second decays in time in accordance 



with Landau damping, consistent with the discussion of 
Linear Landau damping 

For this classical test problem, we choose the usual initial condition fo(x, v) = A cos(kx)fM(v), 
with A = 0.01 and k = 0.5. For the Maxwellian distribution function the dispersion relation 
becomes 

If. oj „ / UJ 



e(k,u) = l + — {l + ^Z. _ 
k 2 \ y/2k \V2k 

where the plasma Z-function is defined as 

Z(z) = -j= r e- t2 ^- z = 2ie- z2 /" e^dt. 



From this relation, the predicted damping rate is computed to be 0.153359 and the predicted 
oscillation frequency to be 1.41566. 

In Figures [4], we plot the evolution of the maximum of the electric field E max using various 
polynomial spaces. In Table [2j we compare the theoretical and numerical values of damping 
rate and frequency as a measurement of accuracy. We see that refining the mesh always 
gives better approximations. The piecewise constant polynomials P° give much larger error 
compared to higher order polynomials. While the difference between the P' and Q l spaces is 
not significant. Observe from Figure [4] how similar the recurrence behavior is for this LVP 
problem to that of the advection equation. Since the linearized equation involves an operator 
Ef' (v), where the electric field depends on the distribution function / on all cells, it is not 
trivial to generalize the analysis of Section|3]to this case. However, it was proven in [381135] that 
there exists a generalization of the Hilbert transform that maps the solution of the advection 
equation to the solution of this LVP system, so this close relation can be justified. 
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Table 2: The damping rate and frequency for linear Landau damping. The numerical values 
are computed using the fourth to the tenth peak and the predicted value is obtained from the 
plasma dispersion function ( 15 ) with a Maxwellian equilibrium. 





Predicted value 


Mesh 


P° 




p2 


Q 1 


Q 2 


Damping rate 


0.153359 


40 x 40 


0.227489 


0.153536 


0.153375 


0.153425 


0.153379 


80 x 80 


0.191702 


0.153366 


0.153363 


0.15369 


0.153363 


Frequency 


1.41566 


40 x 40 


1.38249 


1.41643 


1.41643 


1.41643 


1.41643 


80 x 80 


1.40056 


1.41576 


1.41576 


1.41576 


1.41576 



As for conservation properties, the charge and momentum are well conserved as predicted 
by Propositions 1 and 3. However, the linear energy Hl demonstrates different behaviors 
depending on the polynomial spaces. Figure [5] shows that Hl decays significantly for all P* 
spaces even upon mesh refinement. On the other hand, the seems to conserve it much 
better. We note that Q 1 conserves Hl much better than P 2 , although the former is a subspace 
of the later. 

Also, note from Figure [6] that the electrostatic energy for both choices of polynomial 
spaces damps at a rate given by twice the Landau damping rate. This is to be expected for 
the linear theory, since after integration over space the oscillatory component is removed and 
E ~ exp(— 2^t). Therefore, if the energy is conserved numerically this damped electrostatic 
energy must be converted into the relative kinetic energy that is represented by the first term 



of (16). Thus, conservation of Hl serves as a global measure of the ability of an algorithm to 
resolve fine scales in velocity space. That this transference must take place for the linear VP 
system was proven in Section IV of Ref. |38j . 



4.2 Nonlinear VP system 

In this section, we consider the nonlinear VP system. As noted above, we benchmark the 
solver against three test cases: the nonlinear Landau damping, two-stream instability, and an 
external drive problem with dynamically accessible initial condition. 

The n-th Log Fourier mode for the electric field E(x, t) (23] is defined as 



logFM n {t) = log 10 




E(x, t) sin(knx) dx 



to 




+ 


I 



-L 

E(x, t) cos(knx) dx 



We will use this quantity to plot data from our various runs. 



Nonlinear Landau damping 

For this case we choose fo(x,v) = /m(^)(1 + Acos(kx)) with A = 0.5, k = 0.5, L = 4tt, and 
V c = 6. We implement the scheme on a 100 x 200 mesh and integrate up to T = 100 using 
three methods: P 2 , P 2 with the positivity-preserving limiter, and Q 2 . In Figure [7], we plot 
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the evolution of the first four Log Fourier modes as a function of time. All three methods 
give qualitatively similar results that compare well with other calculations in the literature. 
We observe initial damping (until t ~ 15), followed by exponential growth (until t pa 40), and 
finally saturation of the modes. Note the predicted recurrence times T R for each of the modes 
are as follows: for logFM u T R = 209.44; for logFM 2l T R = 104.72; for logFM 3l T R = 69.81; 
and for logFM^, Tr = 52.36. Since the bounce time is about 40, we have some confidence that 
the solution is resolved at least up u ntil nonlinearity becomes important. Although, the role 
played by T R for the nonlinear evolution is not clear since nonlinearity could remove the fine 
scales generated by linear phase mixing. 



In Figure pi we plot the conserved quantities of Section 2.2 The charge and momentum 



are well conserved for all methods, while the enstrophy has decayed by about 15% at T = 100 
for all three methods. This result agrees with our analysis in Section 2. We remark that 
the limiter has an effect on charge conservation, due to its modification of the solution on 
the boundary. The total energy is conserved much better without the positivity-preserving 
limiter. When we use the limiter, the total energy grows by about 0.3% at T = 100. 



Two-stream instability 

For this case we choose fo(x,v) = fTs(v)(l + Acos(kx)), where frsiv) = ~^ v2e ~ v '' '^ \ A = 
0.05, k = 0.5, L = 4tt, and V c = 6. The mesh size we take is 100 x 200. In Figure |9j we 
plot the evolution of conserved quantities. For this example, charge and momentum are well 
conserved by all methods, so are not plotted. The enstrophy decays by about 4% at T = 100, 
while the total energy is well conserved even with the limiter. The plots of the log Fourier 
modes show an early exponential growth followed by oscillation. Figure 10 provides evidence 
that the system has relaxed into a BGK mode. Here, the relation defined by the ordered pair 
(e = v 2 /2 + $(x, T), f(x, v , if:))) is plotted at various times t. The use of this kind of plot 
as a diagnostic was first reported in [23] for electrostatic VP equations and later in [H] for 
the gravitational VP equations. Here, the evolution clearly indicates convergence to a BGK 
equilibrium. 



Dynamically accessible excitations— KEEN waves 

Motivated by experiments performed for understanding aspects of laser-plasma interaction 
[33] . several authors have considered numerical solution of the VP system with a transitory 
external driving electric field (see, e.g., [TJ [25] ), rather than just specifying an ad hoc initial 
condition for /, as is usually done. Such drive generated initial conditions are examples of those 
proposed and discussed in [36[ [3T], EE] , where they were termed dynamically accessible (DA) 
initial conditions. DA initial conditions are important because they have a Hamiltonian origin 
and preserve phase space constraints. Moreover, since ultimately any perturbation of charged 
particles within the confines of VP theory is in fact caused by an electric field, it is physically 
very natural to consider DA initial conditions. We consider two numerical examples and 
compare our results with those of [U [25], w here the authors observed saturation to nonlinear 
traveling BGK-like states that they termed KEEN waves, standing for kinetic electrostatic 
electron nonlinear waves. We note that the calculations of [1] were duplicated in [22] and 
allied work was given in [1H1 SI] . 
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Specifically, the system is driven by a single prescribed frequency and wavelength, where 
the driven Vlasov equation, 

ft + vf x -(E + E ext )f v = 0, 

is solved. Here, E ext (x,t) = A d {t) sin(kx — out) is the external field, where A d is a temporal 
envelope that is ramped up to a plateau and then ramped down to zero. For out two examples 
we consider the following two ramping profiles: 



A J d {t) = { 



A m sinO/100) 



if < t < 50 



A m if 50 < t < 150 

A m cos ((t - 150)tt/100) if 150 < t < 200 







if 200 < t < T 



(17) 



with A m = 0.052 as used in [2S] and 



.4 



i 



A A d {t) 



m i_|_ e -40(t-io) 



A 1 1 

■"-m \ J- i_|_ e -40(t-lio) 



if < t < 60 



if 60 < t < T 



with A m = 0.4 as used in pQ. In practice, the system is initialized on f(0,x,v) = /m( u ), then 
ramped according to ( Jl7| ) or (18) to prepare the DA initial condition. The system is then 
evolved after E ext is turned off and seen to approach asymptotic states. For both cases the 
computational domain is of size [0, 2ir/k\ x [—8,8], and we take k = 0.26 and u = 0.37. 

Following j2S] with the drive A J d of (17) with A m = 0.4 we obtain for latter times a 
translating BKG-like state, a snapshot of which is depicted in the phase space portrait of 



Fig. 11 This structure moves through the spatial domain giving rise to the central periodic 



electric field signal, E(0,t), depicted in Fig. 12 The period of this signal coincides with the 
propagation speed of the BKG-like state, which in agreement with [25] is about 1.35. Figure 
T3l shows the first four Fourier modes and indicates saturation. 

Next, we increase the drive to compare with results of [T]. With the stronger drive of A d 
with A m = 0.4, the system does not approach a uniformly translating state, but approaches 
a structure with more complicated time dependence as seen in the phase contour plots of 



Fig. 14 These figures are in good agreement with those of [T]. 



The electric field in the middle of the spatial domain, E(0,t), is plotted in Fig. 15, which 
shows more complicated behavior, which surprisingly heretofore has not been plotted. In the 
top part of this figure we see that there is regular periodic behavior at long times and from 
the bottom part of the figure we see that there is period-4 modulation of a basic periodic 
structure similar to that of Fig. 12 Closer examination of phase space plots shows that this 
modulation is cause by the existence of additional smaller BGK-like structures. We note, 
that the existence of multiple BGK-like states is not new; for example, they were seen in the 
simulations of [TB]. Thus, we propose that KEEN waves can be interpreted as the interaction 
of multiple BGK-states, which can also be interpreted as an infinite-dimensional version of 
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Lyapunov-Moser-Weinstein periodic orbits in Hamiltonian systems (see, e.g. [39]). This will 
be the subject of a fu ture publication, so we do not pursue it further here. 

Finally, in Fig. 16 we see from the evolution of log Fourier modes. Prior to t = 10 the 
solution remains roughly at Maxwellian equilibrium. However, at around t = 45 we can 



observe the formation of the KEEN wave, which continues to execute the behavior of Fig. 15 



well after the external field has been turned off at t=60. We see from this figure the effects of 
mesh refinement and the use of different polynomial bases, as indicated in the figure. 



5 Conclusion 

In this paper, we considered the RKDG method for the VP system. We focused on two 
common solution spaces, viz., those with P ; and Q 1 elements. Ignoring boundary contributions, 
the scheme can preserve the charge and momentum, and maintain the total energy up to 
approximation errors when the polynomial order I is taken big enough. However, when the 
positivity-preserving limiter was used, some examples gave relatively large errors for the total 
energy. A rigorous study of numerical recurrence was performed for the Q' elements, and the 
eigenvalues of the amplification matrix were explicitly obtained. DG schemes of higher order 
were shown numerically to give a recurrence time that is close to the classical calculation Tr = 
The qualitative behaviors of the P' and spaces were similar for most computational 
examples, except the linear energy Hl was much better conserved usi ng the Q l space. The 
schemes were used to compute the test cases of Landau damping, the two-stream instability 
and the KEEN wave, and results comparable to those in the literature were obtained. 
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(c) Maxwellian, Q 2 (d) Lorentzian, Q 2 

Figure 2: Computations of the advection equation for the polynomial spaces Q 1 and Q 2 
showing local maxima of the density p ma x as a function of time. The mesh is 40 x 40 with 
Ax = tt/10. For the Maxwellian equilibrium Av = 1/4, while for the Lorentzian equilibrium 
Av = 3/2. 
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(c) Maxwellian, P 2 (d) Lorentzian, P 2 

Figure 3: Computations of the advection equation for the polynomial spaces P 1 and P 2 showing 
local maxima of the density p max as a function of time. The mesh is 40 x 40 with Ax = 7r/10. 
For the Maxwellian equilibrium Av = 1/4, while for the Lorentzian equilibrium Av = 3/2. 
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Figure 4: Depiction of linear Landau damping showing recurrence in the maxima of the electric 
field, E max , as a function of time for various polynomial spaces. 





(a) 40 x 40 mesh 



(b) 80 x 80 mesh 



Figure 5: Evolution of the linear energy Hl of (16) as a function of time, while the Vlasov 



system undergoes linear Landau damping. Various polynomial spaces and mesh sizes were 
used, as indicated. 
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Figure 6: Evolution of the electrostatic energy (red), linear energy (blue) and the first term 
in the linear energy (green) as a function of time, while the Vlasov system undergoes linear 
Landau damping. Here Q 2 was used with a 80 x 80 mesh. 
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Figure 7: Evolution of the first four Log Fourier mode as a function of time for nonlinear 
Landau damping. Various values of the numerical damping/growth rate are marked on the 
graphs. Here the P 2 space with the positivity-preserving limiter was used on a 100 x 200 mesh. 
The predicted recurrence time Tr for logFMi is 209.44, for logFM 2 is 104.72, for logFM 3 is 
69.81, and for logFM 4 is 52.36. 
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Figure 8: Evolution of conserved quantities as a function of time during the course of nonlinear 
Landau damping for various computational methods. A mesh of 100 x 200 was used. 
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Figure 9: Depiction of the first four Log Fourier modes during the nonlinear evolution of the 
two-stream instability. Also depicted is the evolution of energy and enstrophy as a function 
of time for various methods. 
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(c) t = 100 (d) t = 200 



Figure 10: Plots of the distribution f(x,v,t) versus e = ^ — for the two-stream 

instability at the times t indicated, showing saturation to a BGK state. Here the P 2 space 
with the positivity-preserving limiter was used on a 100 x 200 mesh. 
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Figure 11: Phase space contour at T = 1000 for with DA initial condition with drive A d . 
The plot suggests saturation to a moving BGK-like state. Here the Q 1 element was used on 
a 200 x 400 mesh. 




Figure 12: The electric field E(0,t) at the center of the spacial domain at late times for the 
drive A J d . The periodicity matches the propagation of the BJK-like state through the domain. 
Here the Q 1 element was used on a 200 x 400 mesh. 
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Figure 13: The first four Log Fourier Modes for the drive A J d , indicating saturation. Here the 
O 1 element was used on a 200 x 400 mesh. 
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(j) t = 160 (k) t = 225 (1) f = 300 

Figure 14: Phase space contour plots for the KEEN wave at the times indicated. A large 
amplitude drive of A m = 0.4 was used, along^with the P 2 basis and a posit ivity-preserving 
limiter on a 200 x 400 mesh. 
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Figure 15: (top) The electric field E(0,t) at the center of the spacial domain at later times for 
the drive A^. The periodic structure is due to multiple interacting BGK-like states, (bottom) 
Blow up indicating a period-4 modulation of a central hole such as that of Fig. [l2| . The 
simulation was done with P 2 elements with a limiter on a 200 x 400 mesh. 
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Figure 16: Evolution of the first four Log Fourier modes as a function of time for the drive of 



Eq. (18). The simulation used P 2 elements with a limiter on a 200 x 400 mesh. 
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